//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html

#include "DGDiff.h"

// MOOSE includes
#include "MooseVariableFE.h"

#include "libmesh/utility.h"

registerMooseObject("snakeApp", DGDiff);

template <>
InputParameters
validParams<DGDiff>()
{
  InputParameters params = validParams<DGKernel>();
  // See header file for sigma and epsilon
  params.addRequiredParam<Real>("sigma", "sigma");
  params.addRequiredParam<Real>("epsilon", "epsilon");
  params.addParam<MaterialPropertyName>(
      "kappa", "kappa", "Property name of the  heat conductivity material property");
  params.addParam<std::vector<BoundaryName>>(
      "exclude_boundary", "The internal side sets to be excluded from this kernel.");
  return params;
}

DGDiff::DGDiff(const InputParameters & parameters)
  : DGKernel(parameters),
    _epsilon(getParam<Real>("epsilon")),
    _sigma(getParam<Real>("sigma")),
    _kappa(getMaterialProperty<Real>("kappa")),
    _kappa_neighbor(getNeighborMaterialProperty<Real>("kappa"))
{
  // Gather information on broken boundaries
  std::vector<BoundaryName> bnd = isParamValid("exclude_boundary")
                                      ? getParam<std::vector<BoundaryName>>("exclude_boundary")
                                      : std::vector<BoundaryName>(0);
  auto bnd_ids = _mesh.getBoundaryIDs(bnd);

  _excluded_boundaries.insert(bnd_ids.begin(), bnd_ids.end());
}

Real
DGDiff::computeQpResidual(Moose::DGResidualType type)
{
  Real r = 0;
  if (excludeBoundary())return 0;
  
  const unsigned int elem_b_order = _var.order();
  double h_elem =
      _current_elem->volume() / _current_side_elem->volume() * 1. / Utility::pow<2>(elem_b_order);
  
  switch (type)
  {
    case Moose::Element:
      r -= 0.5 * (_kappa[_qp]  * _grad_u[_qp] * _normals[_qp] +
                  _kappa_neighbor[_qp]  * _grad_u_neighbor[_qp] * _normals[_qp]) *
           _test[_i][_qp];
      r += _epsilon * 0.5 * (_u[_qp] - _u_neighbor[_qp]) * _kappa[_qp]  * _grad_test[_i][_qp] *
           _normals[_qp];
      r += _sigma / h_elem * (_u[_qp] - _u_neighbor[_qp]) *_kappa[_qp] * _test[_i][_qp];
      break;

    case Moose::Neighbor:
      r += 0.5 * (_kappa[_qp]  * _grad_u[_qp] * _normals[_qp] +
                  _kappa_neighbor[_qp]  * _grad_u_neighbor[_qp] * _normals[_qp]) *
           _test_neighbor[_i][_qp];
      r += _epsilon * 0.5 * (_u[_qp] - _u_neighbor[_qp]) * _kappa_neighbor[_qp] *
           _grad_test_neighbor[_i][_qp] * _normals[_qp];
      r -= _sigma / h_elem * (_u[_qp] - _u_neighbor[_qp]) *_kappa_neighbor[_qp]*
       _test_neighbor[_i][_qp];
      break;
  }

  return r;
}

Real
DGDiff::computeQpJacobian(Moose::DGJacobianType type)
{
  Real r = 0;
  if (excludeBoundary())return 0;
  const unsigned int elem_b_order = _var.order();
  double h_elem =
      _current_elem->volume() / _current_side_elem->volume() * 1. / Utility::pow<2>(elem_b_order);



  switch (type)
  {
    case Moose::ElementElement:
      r -= 0.5 * _kappa[_qp]  * _grad_phi[_j][_qp] * _normals[_qp] * _test[_i][_qp];
      r += _epsilon * 0.5 * _phi[_j][_qp] * _kappa[_qp]  * _grad_test[_i][_qp] * _normals[_qp];
      r += _sigma / h_elem * _phi[_j][_qp] * _kappa[_qp] *_test[_i][_qp];
      break;

    case Moose::ElementNeighbor:
      r -= 0.5 * _kappa_neighbor[_qp] * _grad_phi_neighbor[_j][_qp] * _normals[_qp] * _test[_i][_qp];
      r += _epsilon * 0.5 * -_phi_neighbor[_j][_qp] * _kappa[_qp] * _grad_test[_i][_qp] *
           _normals[_qp];
      r += _sigma / h_elem * -_phi_neighbor[_j][_qp] *_kappa[_qp]* _test[_i][_qp];
      break;

    case Moose::NeighborElement:
      r += 0.5 * _kappa[_qp] * _grad_phi[_j][_qp] * _normals[_qp] * _test_neighbor[_i][_qp];
      r += _epsilon * 0.5 * _phi[_j][_qp] * _kappa_neighbor[_qp] * _grad_test_neighbor[_i][_qp] *
           _normals[_qp];
      r -= _sigma / h_elem * _phi[_j][_qp] *_kappa_neighbor[_qp]* _test_neighbor[_i][_qp];
      break;

    case Moose::NeighborNeighbor:
      r += 0.5 * _kappa_neighbor[_qp] * _grad_phi_neighbor[_j][_qp] * _normals[_qp] *
           _test_neighbor[_i][_qp];
      r += _epsilon * 0.5 * -_phi_neighbor[_j][_qp] * _kappa_neighbor[_qp] *
           _grad_test_neighbor[_i][_qp] * _normals[_qp];
      r -= _sigma / h_elem * -_phi_neighbor[_j][_qp] *_kappa_neighbor[_qp]*
       _test_neighbor[_i][_qp];
      break;
  }

  return r;
}
bool
DGDiff::excludeBoundary() const
{
  if (_excluded_boundaries.empty())
    return false;

  auto boundary_ids = _mesh.getBoundaryIDs(_current_elem, _current_side);
  for (auto bid : boundary_ids)
    if (_excluded_boundaries.find(bid) != _excluded_boundaries.end())
      return true;

  // make sure we will also break on the neighboring side
  unsigned int neighbor_side = _neighbor_elem->which_neighbor_am_i(_current_elem);
  boundary_ids = _mesh.getBoundaryIDs(_neighbor_elem, neighbor_side);
  for (auto bid : boundary_ids)
    if (_excluded_boundaries.find(bid) != _excluded_boundaries.end())
      return true;

  return false;
}